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We model the dielectric properties of the relaxor ferroelectric perovskite solid solution BaTiOz 
- Bi(Zn 1 / 2 Ti 1 / 2 )03 (BT-BZT) using ab initio methods combined with an Ising-like model. BT- 
BZT has both A-site and B-site disorder, which makes it a challenging material to model. The 
polarization of all possible 2x2x2 supercells of BT-BZT are determined using ab initio methods. 
The cells are then arranged on a lattice in an Ising-like model to determine long-range dielectric 
properties. This enables us to predict the long-range disorder effects of relaxors, while keeping the 
accuracy of ab initio methods. This new model enables us to work with orders of magnitude more 
atoms than other models such as molecular dynamics. We analyze the Monte Carlo data for a single 
lattice configuration using the multiple histogram method, and develop a new modified histogram 
technique to combine data from multiple lattice configurations. Experiments qualitatively agree 
with our models predicted values of dielectric constant and polarization. 



I. INTRODUCTION 

Relaxor ferroelectrics (commonly known as relaxors) 
feature a broad maximum of their dielectric permittivity, 
and a strong frequency dependence of the temperature 
at which that maximum occurs. Relaxors have been the 
subject of considerable research^, due to their funda- 
mental scientific interest, as well as their use in techno- 
logical applications such as capacitors and piezoelectric 
devices^. 

Burns and Dacota first ascribed the diffuse phase 
transition seen in relaxors to the formation of polar 
nano-regions (PNR) below the so-called Burns tem- 
perature (T;,). The PNRs consist of groups of cor- 
related dipoles, which explains the anomalous relax- 
ation times observed in the frequency-dependent dielec- 
tric response. The PNRs generally appear as a re- 
sult of compositional disorder— 7 ~— , which poses a chal- 
lenge to model effectively. There have been many at- 
tempts with varying levels of success. These methods 
include molecular dynamics (MD) simulations^—, den- 
sity functional theory (DFT) calculations^—, spherical 
random-bond-random-field models, random ficld-Ising- 
models22i2£, other Monte Carlo (MC) simulations 2 ^. For 
a good summary of many of these methods applied to 
relaxors, see Burton et. alJ£. While DFT calculations 
are quite accurate, they are limited to small supercells of 
no more than a few dozen atoms. MD calculations are 
able to observe a few thousand atoms before the compu- 
tational demands become too great. MC simulations are 
able to model enough atoms to capture the long-range 
disorder effects needed to describe relaxors, but do not 
have the accuracy of ab initio methods. We aim to create 
a model using both DFT and a MC Ising-like method to 
try and capture the accuracy of DFT, while utilizing the 
long-range scaling of Ising models. 

In this paper, we study the xBaTi03 + (1 — 
x)Bi(Zn 1/2 Ti 1/2 )0 3 solid solution (BT-BZT), which has 
been recently studied and found to exhibit relaxor be- 
havior, with a striking peak value of permittivity 2 ^—. 



Unlike many relaxor ferroelectrics, BT-BZT exhibits dis- 
order in both the A-site and the B-site. Experimentally, 
BT-BZT is found to be slightly tetragonal for composi- 
tions containing mostly BT, with a c/a ratio of approxi- 
mately 1.01 as an upper limits. Experiments also show 
that as the composition of BT-BZT contains less BT, the 
sharp ferroelectric phase transition is lost, and the solid 
solution transitions to relaxor behavio r 25 ' 27 . 



II. MODEL 

We model this relaxor behavior by using a combina- 
tion of techniques including ab initio methods and an 
Ising-like model. The ab initio calculations arc accurate 
at small length scales, but lack the ability to describe 
long-range interactions. The Ising-like model is able de- 
scribe long-range effects, but cannot predict short-range 
interactions. By combining these methods into a single 
model we gain insight into the larger picture by including 
both long and short-range effects. 

Our ab initio calculations are restricted to 2 x 2 x 2 su- 
percells of BT-BZT for symmetry. This allows for compo- 
sitions of x — 0, 0.25, 0.5, 0.75, 1, in the 40 atom supercell. 
Table [J shows the number of distinct configurations for 
each composition, and Figure[T]shows one possible atomic 
configuration for a 2 x 2 x 2 supercell with a single zinc 
atom [x = 0.75). With the Quantum Espresso^ pack- 
age we use Density Functional Theory (DFT) and the 
Modern Theory of Polarization (MTP)22r^i to determine 
the energies and polarizations of each unique configura- 
tion. The DFT calculations were run using a PBE-GGA 
exchange-correlation functional, ultrasoft pseudopoten- 
tials, a cutoff energy of 80 Rydberg, and a 2 x 2 x 2 
Monkhorst-Pack k-point grid. 

The polarizations determined by MTP are not unique 
for each cell. The actual polarization is only accurate 
modulo eR/Q, where e is the charge of an electron, R is 
a lattice vector, and fl is the cell volume. We are able 
to determine the polarization of the relaxed cell relative 
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to the polarization of a cell with atoms in the symmetric 
positions. Some configurations may also have degenerate 
polarization states. For cells with more than one polar- 
ization, all accessible states are recorded and used later 
in the Monte Carlo simulations. 

There are some disadvantages to using ab initio meth- 
ods to model relaxors. Most importantly, relaxors display 
long-range disorder, which is not efficiently modeled us- 
ing ab initio techniques. By using a lattice model, we are 
able to observe some long-range effects, as well as study 
any compositional value x. 
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FIG. 1: (color online) BT-BZT cubic 2x2x2 supercell ex- 
ample with x = 0.75. The lattice constant of the supercell is 
twice that of the smallest cubic BT unit cell. 



TABLE I: The number of unique atomic configurations for 
each composition of a 2 x 2 x 2 supercell. 



Number of Zn atoms 


Composition x 


Distinct configurations 





1.0 


1 


1 


0.75 


3 


2 


0.50 


26 


3 


0.25 


13 


4 


0.0 
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Once we have determined the energies and all polariza- 
tion states of each configuration for all compositions, we 
are able use this data as the basis of an Ising-likc model. 
The classic Ising model is used to determine critical be- 
havior for ferromagnetic systems on a lattice, where each 
lattice may have a spin up/down state^r— . An exten- 
sion of the Ising model is the Heisenberg model which 
may have spins in any direction. Here, we create a lat- 
tice of our 2x2x2 supercells and use only the polariza- 
tions determined by the MTP calculations. This differs 
from both Ising and Heisenberg models in that the po- 
larizations are not simple up and down, but neither are 
they allowed to point in any direction. This adaptation 
is not unlike a Potts model^— in that there are discrete 
directions of the polarizations, but our polarizations are 



not uniform in nature and each cell on the lattice may 
have different accessible polarization states. We use the 
interaction Hamiltonian of the Heisenberg model 



H = -- 



J 



cells NN 



cells 



E 



(1) 



where Pi is the polarization of a single cell, E is the exter- 
nal electric field and J is a coupling constant. The cou- 
pling constant is tuned to match the experimental Curie 
temperature of pure barium titanate using the fourth or- 
der Binder cumulant method^. 

Using the energies and polarizations from the ab initio 
calculations, we populate a lattice of cells to obtain the 
desired composition x. The number of each unique 2 x 2 x 
2 cells is determined using Boltzmann statistics, and the 
cells are distributed on the lattice stochastically. Once 
the lattice is established, we begin with the Monte Carlo 
simulation. 

For each desired temperature and composition, we 
equilibrate for at least 10 6 flips per cell. The total po- 
larization and energy are recorded after every 10 changes 
per site. Equilibrium is considered to be reached when 
the system has run for several energy correlation times. 
Once the system has reached equilibrium, we can calcu- 
late the mean value of any function of the energy and 
polarization simply by 



1 N 



(2) 



where the subscript i is the index of a particular data 
sample, and N is the total number of samples. We can 
then calculate properties such as the dielectric constant 
and specific heat using the variance of the polarization 
and energy respectively. 

Unfortunately, we need to run a separate simulation 
for each temperature and composition. By using the his- 
togram method, a single simulation can be used to de- 
termine results near the simulation temperature^. This 
idea is expanded to use simulation data at multiple tem- 
peratures to determine results for a large range of tem- 
peratures. 

The multiple histogram equations are given by 



JV 



l g m exp {-P m E - f m ) 



(3) 



where j3 m — 1/ksTm, N is the number of simulations, 
H n (E,p) is the histogram of the n th simulation, n m is the 
total number of samples included in the m th histogram, 
g n = 1 + 2t„ where r n is the correlation time, and f m is 
an estimate of the dimensionless free energy of the m th 
simulation. The correlation time is incorporated as a 
weighting factor to account for the differing correlation 
times of the samples for different simulations. The free 
energies serve as a correction factor due to combining 
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data from multiple temperatures, and are calculated self 
consistently using 



exp(/ m ) = Pp m {E,p) 



(4) 



Note that Equation [3] does not give the actual proba- 
bility, but a probability distribution function that must 
be normalized: 



Pp{E,p) 



(5) 



Once the probabilities are determined, we are able to 
calculate the mean value of a quantity using 



N K N„ 



{f(E,p)) p = E E fiEuPjWEuPi) , (6) 

i 3 

By calculating the mean of the polarization in a given 
direction (P x ), as well as (P x 2 ), we can use the variance 
of the polarization to calculate the susceptibility by 



Xx* = P {(P x 2 ) - {P x ?) 



(7) 



and the dielectric constant by e = 1 + 47rx- 

Figure [2] compares the dielectric constant for single 
simulations at every 10 K and the multiple histogram 
method using the data from the same input simulations. 
The multiple histogram method can be thought of as an 
interpolation scheme if there are enough data samples at 
all of the simulation temperatures. 



the cells on the lattice. If we wish to combine the results 
of multiple arrangements into a single result, we must use 
an appropriately weighted average. The weighted result 
must be equivalent to using one large lattice that con- 
tains each other lattice as a subset. In order to combine 
these results we multiply the probabilities and sum over 
all states with an appropriate Kroncckcr delta. 
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X Pp,N {En,Pn) X 5e,(E 1 + E 2 + ...+E n )' > < 

Sp 

,(Pi+P2 + ---+Pn) ' 

(8) 

where the subscript /3 indicates the dependence of proba- 
bility on temperature. The total probability is a function 
of the sum of all energies and polarizations. This works 
on the assumption that the systems are non-interacting 
and independent. This method is quite expensive com- 
putationally, so we utilized convolutions to speed up the 
calculations. This technique may also be applied to mul- 
tiple runs of the same configuration, to avoid having to 
repeat costly calculations. We are also able to use this 
method to combine polarization data from all cartesian 
directions into a single weighted value. Figure [3] shows a 
single cartesian component of the dielectric constant for 
five different lattice configurations of a single composi- 
tion, along with the weighted combined result. 
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FIG. 2: (color online) Dielectric constant of BT (x — 1) on 
a 16 x 16 x 16 lattice using raw data and the multiple his- 
togram method. The simulations used for the old method are 
also the inputs in this multiple histogram method. The mul- 
tiple histogram method is sampled every IK, while the input 
simulations were performed every 10K. 



A. Combined-Multiple Histogram Method 

Since not every cell on the lattice is identical for com- 
positions of x < 1, there is more than one way to arrange 



FIG. 3: (color online) The dielectric constant of a single com- 
ponent of the polarization for each of five different lattice 
configurations and the combined result for x — 0.95. 

We are able to determine the uncertainty in our cal- 
culations with the regular multiple histogram method by 
first finding the uncertainty in the probabilities: 

6P P (E,p) = (E 1 P P (E,p) , (9) 

Vn=l 9n J 

where the sum is over all simulations, H n (E,p) is the 
histogram of a single simulation, and g n ~ 1 + 2r„ is the 
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correlation time factor—. The uncertainty in the total 
probability of Equation [5] is then given by 



N 



N 



(10) 



where the uncertainty for each lattice configuration SPi 
is calculated using Equation [3] 

The derivative of Equation [8] with respect to a single 
lattice configuration probability is 
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(12) 



The uncertainties of the susceptibility can be deter- 
mined using 



*x 2 = £ 



i\dPp(E,p) 



6P p (E,p) 



(13) 



By carrying out the derivative, we get 

*x 2 = ^ E [t^ 2 - <p 2 > + 2 ^ 2 - 2 <*M ! 



where 



{5Pp{E,p)f , (14) 



Z p = J2Pp(E,P). (15) 



III. RESULTS AND DISCUSSION 

A. Dielectric Constant 

The results of our model are compared with the ex- 
periments of Huang and Cann2£, where all of the results 
produced by the model were found using the combined 



histogram method. Each combined data set is found us- 
ing simulations performed on a 16 x 16 x 16 lattice at 
every 5 Kelvin, each with at least 10 6 Monte Carlo flips 
per cell. 

Figure 0] shows the dielectric constant for the model 
and experiment for various compositions. Each of the 
experimental results of the dielectric constant are shown 
for a range of frequencies. In general, at high tempera- 
tures there is little to no frequency dependence of the ex- 
perimental dielectric constant, while at low temperature 
there is some frequency dependence that grows larger as 
the composition decreases. The predicted dielectric con- 
stant at low temperatures does not match experiment for 
any composition, due to the frequency dependence and 
the fact that our model only predicts the zero-frequency 
limit. 

The top section of Figure 2] shows the dielectric con- 
stant for pure barium titanatc (BT) (x — 1). The experi- 
mental results show a clear transition near the Curie tem- 
perature of 393-ftT, and essentially no difference between 
the various frequencies. The model prediction shows a 
very precise curve, as the error bars are smaller than the 
thickness of the line. However, the value of the curve is 
not as accurate as one could hope. The predicted dielec- 
tric constant is well above the experimental value for all 
temperatures. For temperatures slightly above the crit- 
ical point, the dielectric constant is almost an order of 
magnitude too high, and at high temperatures the differ- 
ence is around a factor of two. 

Below the critical point, the predicted value diverges, 
while the experimental value drops down significantly. 
This divergence can be understood by considering that 
the model computes the zero frequency limit (DC) in 
thermal equilibrium^!. 

At low temperatures, a ferroelectric is not likely to 
change polarizations over short time scales. Thus, only 
a single polarization state is all that is needed to deter- 
mine the dielectric properties. The combined histogram 
method includes many polarization states for each tem- 
perature, and thus gives a non-physical result below the 
critical temperature. One possible solution would be to 
break the symmetry of the system by calculating the 
probabilities used to determine the dielectric constant in 
the presence of a small applied electric field. This would 
force the system to have only those polar states that are 
aligned with the electric field. In practice, the symme- 
try breaking leads to instabilities in both the dielectric 
constant and specific heat, which is likely due to poor 
convergence in our low temperature simulations. 

The middle section of Figures 0] shows the experimen- 
tal and model data for BaTi0 3 - Bi(Zn 1/2 Ti 1/2 )0 3 (BT- 
BZT) with a composition of x = 0.95. The model data 
was found using the combined histogram method with 
the x, y, and z components of the polarization for five 
different lattice configurations. This is roughly equiva- 
lent to using a single component of the polarization for 
15 separate lattice configurations. Since there are mul- 
tiple configurations being considered, there is a visibly 
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FIG. 4: (color online) The dielectric constant of three com- 
positions, x = 1, x = 0.95 and x = 0.8. The predicted result 
is compared with experimental measurements made at sev- 
eral frequencies^. Theory is in black, with error bars given. 
Experiment is shown with colored symbols. 



larger uncertainty in the dielectric constant for x = 0.95 
than for pure BT (x = 1). 

The experimental data for x = 0.95 still shows some 
critical behavior at a slightly lower temperature than 
pure BT. The model predicts the dielectric constant rea- 



sonably well at high temperatures. The difference be- 
tween the model data and experimental data is now much 
less than a factor of two, even as the temperature de- 
creases. For lower temperatures the model predictions 
still diverge for these compositions. The temperature 
at which this divergence occurs, however, is much lower 
than that of pure BT. 

The lower section of Figures 3] shows the experimen- 
tal and model data for BaTi0 3 - Bi(Zn 1/2 Ti 1/2 )0 3 (BT- 
BZT) with a composition of x — 0.8. These model re- 
sults were found using the three polarization components 
of only three atomic configurations. The model dielectric 
constant for this composition is slightly lower than exper- 
iment for nearly all temperatures, while the actual values 
are very similar. 

B. Non-Zero External Electric Field 

So far our discussion and analysis has only included 
effects from the nearest neighbor interaction term in the 
Hamiltonian in Equation [TJ For ferroelectric materials 
it helps to understand the behavior of the material in 
the presence of an external electric field. To include the 
second term in the Hamiltonian in our analysis, we only 
need to modify our probability (Eq. [3]) to include this 
new term in the energy: 

E»=i H n {E,p)g- 1 exp (-/3 (e -p • e)) 
E m =i n ™9™ exp (-p m E - j m ) 

(16) 

Note that the energy in the denominator does not in- 
clude the electric field term, because all simulations were 
done with zero electric field applied. The beauty of the 
histogram method is that as long as the polarization and 
energy microstatcs exist in the histogram, then the prob- 
ability of finding the system with that microstate can be 
found. The energy in the denominator of Equation [16] 
corresponds to the Hamiltonian used in each simulation. 
This is made clear since it is this energy that is multi- 
plied by /3 m , which corresponds to the simulation tem- 
peratures. Since the value of (3 in the numerator corre- 
sponds to the desired temperature, the desired electric 
field term also makes an appearance. Thus, using the 
histogram method it is possible to use data simulated at 
zero applied E-ficld to get a probability at a non-zero 
E-field. 

The top section of Figure [5] shows the polarization den- 
sity in the direction of the field as a function of electric 
field for both the model and experimental data of pure 
BT at room temperature. The obvious difference be- 
tween the model and experiment is that the model is com- 
pletely polarized in the direction of the applied field, even 
for small fields. Since pure BT is ferroelectric at room 
temperature, the large polarization is expected, even in 
the absence of an applied electric field. 

The magnitude differs from the experimental values at 
low fields because the Monte Carlo model predicts a DC 



6 




-15 

-20 | 

-40 -30 -20 -10 10 20 30 40 
Efield ( kV/cm ) 



FIG. 5: (color online) Polarization vs electric field for various 
compositions compared with experiment at room tempera- 
ture. Experimental data provided by Huang and Canr*2£. 

value, which corresponds to the infinite time limit of an 
applied DC field. This leads to the drawback that we 
are unable to predict the hysteretic behavior shown by 
experiment. The slope of the polarization versus elec- 
tric field curve is the susceptibility in the limit of small 
electric fields. The discontinuity of the predicted curve 
shows an infinite slope for small fields. This corresponds 



with the divergence of the dielectric constant below the 
Curie temperature of BT. In the limit of high field, the 
experiment and model agree quite nicely. This is due to 
the fact that high field strengths force the system into 
a single bulk polarization state that is equivalent to the 
equilibrium state for the model. 

The middle and bottom portions of Figure [5] show the 
polarization density vs electric field with compositions of 
0.95 and 0.8 respectively for both experiment and model 
data. As was previously stated, and may be more obvious 
here, the model fails to predict hysteresis, since there can 
only be one true equilibrium state in the presence of an 
electric field. 

One feature that occurs for both compositions of a; < 1 , 
is that the maximum polarization of the model is much 
lower than the maximum for the experiment. This is 
most likely due to the lack of available microstates with 
high polarizations in the model. Since the input simu- 
lations are all performed with zero applied electric field, 
the high polarization states most likely have not been 
sampled. If the solid solutions had high polarization 
states available as input to the histogram method, the 
maximum polarization of the P vs E graphs would bet- 
ter match the experimental values, which suggests that 
we would benefit from performing additional MC simu- 
lations with an applied finite electric field. Even so, the 
combined multiple histogram method helps us by reduc- 
ing the number of such simulations required. For BT this 
issue did not arise because the high polarizations were 
sampled, since the low temperature states with zero field 
are identical to states with a large electric field applied. 

For all of the solid solution polarization versus electric 
field graphs there is a finite slope at small applied field, 
while the slope of the x — 1 graph is infinite. The lin- 
ear behavior of this slope at small applied fields indicates 
that the materials with lower composition are linear di- 
electrics, and the value of the slope and the dielectric 
constant show that they are in fact high permittivity lin- 
ear dielectrics. 



IV. CONCLUSION AND SUMMARY 

Using only ab initio methods and the experimental 
Curie temperature of BaTi03 (BT), we developed a 
model to predict the dielectric properties of BaTi03 - 
Bi(Zn 1/2 Ti 1/2 )0 3 (BT-BZT) in solid solution. This ma- 
terial behaves as a relaxor ferroelectric for some compo- 
sitions as shown experimentally^ and through our own 
predictions. 

We first use Density Functional Theory (DFT) to cal- 
culate the ground state energies of small 2x2x2 cu- 
bic supercells of the BT-BZT solid solution. This is a 
slight challenge, since there are many atomic configura- 
tions of the solid solution supercells. We then use the 
Modern Theory of Polarization (MTP) to determine the 
available polarization states of each supercell. DFT and 
MTP alone are unable to effectively model relaxors due 
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to the importance of long range disorder. 

In order to model the long-range disorder effects of 
rclaxors we introduce a Monte Carlo lattice Ising-likc 
model, where the cells on the lattice are the supercells 
used in the DFT calculations with polarization states de- 
termined by MTP. Our model differs from Ising, Potts, 
and Hciscnbcrg models in that each cell has a finite num- 
ber of polarization states that arc of differing magnitude 
and orientation, and not every cell on the lattice has the 
same available states. The interaction energy coupling 
constant needed for this Ising-like model was determined 
by fitting to the experimental Curie temperature of BT. 

To analyze the Monte Carlo data we use the multiple 
histogram method to combine data from simulations of 
a single lattice configuration at multiple temperatures. 



However, to combine simulations of multiple lattice con- 
figurations we had to develop a modified version of the 
multiple histogram method. Using this new analysis tool, 
we combine Monte Carlo data from multiple simulations 
and temperatures to obtain the most statistically accu- 
rate data possible. 

The new model and analysis method allow us to pre- 
dict a dielectric constant that has reasonable agreement 
with experiment, especially at high temperatures. The 
manner in which the predicted dielectric constant as a 
function of temperature changes with respect to the com- 
position is consistent with that of experiment. We also 
predict the behavior of the polarization of BT-BZT in 
the presence of an applied external electric field, with 
reasonable agreement in the range of small electric fields. 
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